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Abstract. A challenging difficulty in solving the radial Dirac eigenvalue problem numerically 
is the presence of spurious (unphysical) eigenvalues among the correct ones that are neither 
related to mathematical interpretations nor to physical explanations. Many attempts have been 
made and several numerical methods have been applied to solve the problem using finite element 
method (FEM), finite difference method (EDM), or other numerical schemes. Unfortunately 
most of these attempts failed to overcome the difficulty. As a FEM approach, this work can be 
regarded as a first promising scheme to solve the spuriousity problem completely. Our approach 
is based on an appropriate choice of trial and test functional spaces. We develop a Streamline 
Upwind Petrov-Galerkin method (SUPG) to the equation and derive an explicit stability pa- 
rameter. 



Introduction. 

Studying the properties of electrons in atoms is governed by the Dirac equation which gives 
a complete picture of the electron behavior by means of specifying its energies (eigenvalues) 
in orbitals around the nucleus. Up to date, computing the eigenvalues of an electron in the 
many-electron system (in some methods) is based on determining the corresponding eigenvalues 
in a single-electron system (Hydrogen- like ions), where the eigencouples are used as a basis to 
approximate the electron energies in the entire system. Unfortunately, computing the eigenvalues 
of the electron in the Hydrogen-like ions by numerical methods is upset by the presence of 
spurious solutions (eigenvalues do not match what is physically observed) . The spurious solutions 
annoy the computations, they disturb the solution in a way it becomes no longer reliable. At the 
time, one can identify the spurious eigenvalues, but there is no efficient method to just remove 
them from the entire spectrum without affecting the genuine values. 

The presence of spurious values in the spectrum of the radial Dirac equation and other prob- 
lems has been addressed in most of numerical computations. In [15], the occurrence of the 
spurious roots has been related to incorrect balancing of large and small components spaces H-^ 
and H^', and has been restricted to the positive quantum number k. In solving Dirac equation 
by a mapped Fourier grid [l], spurious values have been detected for k > 1, where the causality 
is recounted to the symmetric treatment of the large and small components. For eigenvalue 
problems in general [2T], the occurrence of spectrum pollution has been related to the absence 
of suitable constraints in the mathematical formulations or discretization, which results in mis- 
matching of desired physical properties of the problem. Shabaev and Tupitsyn et al. [17^ [T9] 
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have also allied the presence of spectrum pollution to the symmetric discretization of the small 
and large components of the wave function. They have pointed out that using the same finite 
space for both components is the essence of the problem. They have proposed an alternative 
method to handle the difficulty by an addition of suitable terms to the basis functions known as 
a basis correction. Also they have explained the property of energy coincidence for the positive 
and the corresponding negative values of k. 

The spuriousity of the eigenvalues computation using spectral Tau method has been studied 
in [6]. Also the causality of the spurious solution in the electromagnetic problems in general 
can be found in [13] . To the Dirac eigenvalue problem, we refer respectively to [20] and [U [T7] 
for finite difference and B-splines approximations. For a brief Finite Element derivation for the 
Dirac operator see e.g [Ti] . 

In the present work, we provide a finite element scheme for solving the radial Coulomb-Dirac 
operator that provides a complete treatment of the spurious eigenvalues. This scheme may be 
considered as the first stable finite element approach for the numerical approximation of the 
Dirac eigenvalue problem. To proceed, we relate the occurrence of spuriousity to the function 
spaces in the implemented numerical method. What ever the method is, finite element method 
(FEM), finite difference method (EDM), the spectral domain approach (SDA), the boundary 
element method (BEM), or the point matching method (PMM), the spuriousity persists. Hence, 
it is priorly understood as not an effect of the numerical method applied, but to a mismatching 
of some physical properties of the eigenstates in the computations. The present work interprets 
the existence of spurious values and their remedy by means of the following two steps: 

(1) The choice of suitable trial functional space that meets the physical property of the wave 
functions in the implemented numerical methods and its role of spuriousity elimination. 

(2) The choice of weighted test functional space, this treats what remains of spurious val- 
ues in one hand, and solves the coincidence of energies for positive and corresponding 
negative k in the other. 

In other words, we classify the spurious solution of Dirac eigenvalue problem into two categories. 
The first is those that appear within the spectrum for all values of k. We call this type the 
instilled spurious values. It is worth to mention that this type of spuriousity appears not only 
for positive k, but for negative k as well. Instilled spurious values affect the true values or may 
degenerate with them which results in some perturbed eigenfunctions. However, this will be 
discussed in detail in the coming section, where, by means of choosing an appropriate space 
of discretization, part of the instilled spurious values is treated. The second category can be 
understood as the coincidence of the first eigenvalue of the radial operator with positive k to that 
with the corresponding negative k. We call this type of spuriousity the unphysical coincidence 
phenomenon: The eigenvalues with positive k have been shown in the finite dimensional spaces 
to be a repetition for those with the corresponding negative k [19], which is not the case in the 
usual (infinite) space of the wave functions. In an attempt to overcome the difficulty, the last 
(main) section is devoted to set a scheme that removes the spuriousity for both categories. 
For a brief sketch of the scheme, consider the radial Coulomb-Dirac equation 
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Here m and c are respectively the mass of the electron and the speed of light, the quantum 



are the total and the orbital angular momentum numbers respectively, and Dx is the derivative 
with respect to x in M. The Coulomb potential, V{x), is given by — , where Z £ [1 , 137] is 
the electric charge number. The unknown / and g are the large and small components of the 
eigenfunction ip with corresponding eigenvalue A. 

The presence of convection terms in the off diagonal and the absence of diffusion terms cause 
numerical instability while computing the eigenvalues A. Indeed, in the standard Galerkin 
finite element solution of the equation one encounters spurious eigenvalues. In order to re- 
move the spuriousity, we derive a stable finite element scheme based on appropriate choice 
of functional spaces of the Dirac spinors. By rewriting the explicit equations of / and g 
and applying suitable boundary conditions, the original space of the Dirac wave functions is 
'Ko{Q,) = {v G C^{U) n Hq{Q) : v'\qq = 0}, where is the space of continuous functions 
which possess continuous first derivatives, Q is some open bounded domain, and Hq{Q,) = {v : 
V and v' are elements of L^(r2), and v\qq = 0} (for all values of k except ±1, where for k = ±1 
the lower boundary condition of v' should differ from zero, but for generality and for sake of 
simplicity it is assumed to vanish, see Remark 1 below). Thus, by this definition, IKo(r2) is the 
space of continuous functions, v, which admit continuous first derivatives that are vanishing 
smoothly on the boundaries. 

Consider the weak form of the equation above of finding {A, (^} G M x !Kq{Q,)'^ such that 



Where D is a test function, and the superscript t is the usual matrix transpose. Cubic Hermit- 
ian (CH) interpolation functions turn out to be a suitable choice which sufficiently fulfill the 
requirements of !Kq{Q,). Let be the finite dimensional subspace of 5{o(il) on the partition kh 
spanned by the piecewise CH basis functions. Choosing t) € (V^)^ as {v,Oy and {0,vY, where v 
is an element of , and assuming f ,g £ , remove partially the first category of spuriousity 
(only for very small Z) and do not help in solving the coincidence phenomenon. 

A complete treatment is achieved by letting the test function to live in another space different 
from that of the trial function, mainly by assuming D to be {v,tv'Y and {tv',vY (where v' 
means D^v) respectively in the variational form above. The scheme is accomplished by deriv- 
ing the stability parameter r, which turns out to have the form r := tj = ^hj^i^j^^^^^^Y 
The derivation is based on two leading simplifications; to consider the limit operator in the 
vicinity of x at infinity (i.e to consider the most numerically unstable part of the operator) and 
c-correspondence dominant parts of the system. From the weak form with the modified test 
function, and after applying the above simplifications we obtain an approximation A(t) of the 
accumulation eigenvalue. Knowing that the limit point eigenvalue is mc^, we like to minimize 
the error |A(t) — mc^|, which gives the desired formula of tj. 

As a numerical method implemented in this work, the finite element method (FEM) is applied, 
with the usual continuous Galerkin method in the first section and a Petrov-Galerkin method in 
the second section. For the integrals evaluation, four-point Gaussian quadrature rule is applied. 



number k is the spin-orbit coupling parameter defined as k 



(_iy+'+2(j + 1)^ where j and / 




The paper is arranged as follows; In the first part we discuss the first category of the spurious 
values and how to remove them partially via choosing suitable trial functional space. A com- 
parison is performed between the incorrect and the correct functional spaces through numerical 
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examples. In the second, we discuss the completion of the treatment. Basically we impose the 
weighted test function to live in a space different from that of the trial function. This is the 
well-known Streamline Upwind Petrov-Galerkin (SUPG) method [21 [71 [10]. Finally a stability 
parameter is derived to achieve the desired goal. 

1. Trial functional space 

Recall the radial Dirac equation 

( mc^ + Vix) c(-D, + f) 
^ ' V c(l)x. + f) -mc^ + V{x) 

where, then, the two-equation system is 

(2) (mc2 + V{x))fix) + c( - g'ix) + ^gix)) = \f{x) , 

(3) c{f{x) + ^/(x)) + ( - + V{x))g{x) = \g{x) . 

We first solve this system by usual continuous linear basis functions (hat functions). Since x 
ranges over [0, cxj), x = represents a singularity for the Coulomb potential and hence careful 
treatment is needed, i.e one can consider extended nucleus in the entire domain or just assume 
point nucleus on a cut-off domain. However, computationally, almost the same technique is used 
for both cases. For simplicity we will treat point nucleus model in all computations except in 
the last table where we extend the computations to extended nucleus. 

Divide the domain Q = [a , b] into n + 1 subintervals with n interior points distributed 
exponentially, and assume : a = xq < xi < X2 < • • • < Xn+i = b the resultant partition of 
with mesh size hi = Xi — 

The exponential distribution of the nodal points is crucial for solving the radial Dirac equation 
in order to get more nodal points near the singularity (x = 0). This is because the wave function 
oscillates much more near the nucleus which means more information is needed about its behavior 
near that region, whereas the fine grid is not required in a position away from the nucleus. 

The choice of the computational space V is important and plays the most influential role in 
the core of the problem. To see that, let us first take the space of only continuous functions as 
the functional space V. We will show, by means of numerical examples, how this space results 
in the occurrence of spurious values. The presence of spuriousity is due to the fact that the only 
continuous functional space lacks to a certain constraint in the mathematical formulation, i.e it 
fails to have an identified property which being exist for the original wave function. 

For a fast and simple algorithm, continuous linear basis functions are considered. So let V = V' 
be the subspace of continuous linear polynomials (the superscript / denotes for the linear case) , 
and let Vj^ C be the finite subspace consists of piecewise continuous linear polynomials on kh 
spanned by the usual linear functions. We assume that both trial and weighted test functions 
belong to this space. For f{x) and g{x) in we write 

n 

(4) f{x) = Y,CM^), 

n 

(5) g{x) = ^ij(l)j{x) , 



j I g{x) -\ g{x) j 
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where Cj and are the unknown values of the functions / and g at the nodal point Xj respectively, 
and is the basis function. Since the eigenfunction decays in the vicinity of x at infinity and 
also considered to be zero at x = 0, the Dirichlet conditions are assumed to treat the boundaries. 
The problem is now read as solving ([2]) and ([3]) such that / = and g = at x = a,b (i.e 
Co = Cn+i = '^0 = ^n+i = 0). The usual finite element method of the problem is to assume / 
and g as above in ([2]) and ([3]), then multiply by a test function and integrate over the domain O 

n n n 

(6) J2 {w^{x)4>j{x) , v{x))Cj + ^{- c(l)'j{x) + -4>j{x) , v{x))^j = {(t>j{x) , v{x))Cj 
and 

n n n 

{w {x)cl)j{x) , v{x))^j = {<t>jix) , v{x))(,j , 
j=i j=i j=i 

where w^{x) = ±mc^ + V{x), and {u , = f^uvdx is the usual L'^{Q) inner product. The 
basis function (pj{x) has its support in [xj^i , Xj] =: Ij and [xj , Xj^i] = Ij+i and defined as 



<f)j{x) 



X G Ij+i 



Let V = he an element of the same space in ([6]) and ([7]), this leads to the symmetric 
generalized eigenvalue problem 

(8) AX = \BX . 

Here A and B are both symmetric block matrices defined by 



(9) A 
and 



mc^Mooo + M( 



V 
000 



cMoio + ckMooi 



-cMoio + ckMooi 



-mc^Mooo + 



000 



(10) B 



where M^^j are n x n matrices defined as 



' Mooo 








Mooo 



(11) {M^,,),, = 4^ x~' q{x) dx , [4\x) = £^Hx)) . 

The vector X is the unknown defined as (C , ^)*, where 

C = (Cl;C2,--- >Cn) 

and 

C = (6,6, ■ ■ ■ ,6) • 

Clearly the diagonal matrices of A are symmetric and the off diagonal matrices consist of two 
parts, one is symmetric and exists in both sides, and the other, (Moio)* = — Mqio, is anti- 
symmetric and exists in both off diagonal sides with different sign, this explains the symmetry 
of the block matrix A. For the block matrix B the symmetry is obvious. 

In Table 1 the first six computed eigenvalues for the Hydrogen atom (Z = 1) are listed for 
\k\ = 1, these eigenvalues are obtained using n = 100 interior nodal points. The exact solution 
for K = — 1 is shown in the right column of the table. Even with mesh refinement the spuriousity 
is still present, see Table 2 with n = 400. 
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Table 1. The first six computed eigenvalues for the electron in the Hydrogen 
atom using linear basis functions with 100 nodal points. 



Level 


K = 1 


K = -1 


Rel. Form, k = —\ 


1 


-0.50000665661 


-0.50000665659 


-0.50000665659 


2 


-0.12500414297 


-0.12500414298 


-0.12500208018 


3 


-0.05556140476 


-0.05556140479 


-0.05555629517 


^ 


^ -0.03192157994 


-0.03192157993^ 


Spurious Eigenvalue 


4 


-0.03124489833 


-0.03124489832 


-0.03125033803 


5 


-0.01981075633 


-0.19810756319 


-0.02000018105 



Table 2. The first six computed eigenvalues for the electron in the Hydrogen 
atom using linear basis functions with 400 nodal points. 



Level 


K = 1 


K = -l 


Rel. Form, k = —1 


1 


^.50000665661 


-0.50000665659 


-0.50000665659 


2 


-0.12500208841 


-0.12500208839 


-0.12500208018 


3 


-0.05555631532 


-0.05555631532 


-0.05555629517 


^ 


gr03l4ll72061 


-0.03141172060"^ Spurious Eigenvalue 


4 


-0.03118772526 


-0.03118772524 


-0.03125033803 


5 


-0.01974434510 


-0.01974434508 


-0.02000018105 



In the tables above, the shaded left corner value is what meant by the unphysical coincidence 
phenomenon, and the values in the fourth row are the so-called instilled spuriousity. The spurious 
values appear for both positive and negative values of quantum number k, and they persist 
despite of mesh refinement. Generally this kind of spurious solution can be identified among 
the right spectrum, but there is no way to just exclude them as a hope of treatment, since they 
have already affected or degenerated the true values. 

As we mentioned before, the occurrence of spuriousity is related to the implementation of 
the numerical method, where the numerical scheme we assumed is the FEM with the proposed 
space VK Therefore, either of them holds the responsibility of causing the spectrum pollution. 
At this end, it is worthy to mention that other methods like finite difference method (FDM), the 
method of moments (MoM) |151 [16] and others, reported the occurrence of spuriousity in many 
computations for the Dirac operator or else. So we conclude that the problem of spuriousity 
is almost caused by the finite element spaces employed in the discretization, and hence the 
causality of spuriousity is V'-problem and never FEM-problem. 

We return to ([2]) and ([3]), rewrite both equations to obtain explicit formulae for / and g 

(12) w+{x){w-{x) - \ff{x) - —{w-{x) - A) (c/'(x) + —fix)) + c\{w~{x) - A) x 



(cfix) + -fix) 

X 



CK. 



X^ 



fix)) - V'ix)icf'ix) + -fix)) = \iw-ix) - Xffix) 



and 
(13) 



w (x) (t(;"'"(rE) — A)^(7(x) -I ( 

X 

icg"ix) - -g'ix) + ^g(x)) - V'ix)icg'ix) - -g(x)) = A(u;+(x) - A)^g(x) . 



[W [X 



\)icg'ix 
^)[cg'i. 



—gix)) + c 

X ' 



ivS^ix) — A) X 



CK 



X 
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Equations ()12p and ()13p can be written in simpler forms as 



(14) fix) + 71 (x, X)f'ix) + 72(x, A)/(x) = 



(15) g"{x) + Oiix, \)g'{x) + 02(x, A)5(x) = . 



Where 



^ {w+{x)-\){w {x)-\) k'^ + k 
12{X, A) = '-^ ^ 




x{w (x) — a) 



kV'{x) 




and 




(w+(2;) — A) [w (x) — A) K? — K 

g2 



3;(fi;+(x) — A) 



kV'{x) 



The terms /" and g" in (jl4p and (jl5p propose further constraint on both components of the 
wave function. By these equations / and g are imposed to be twice differentiable. This means 
that / and g should be continuous with continuous first derivatives, hence the proposed original 
domain is C\n)nH^{n). 

Instead of regarding V' as the space of variation, a space of continuous functions with con- 
tinuous first derivative is considered to discretize both components of the wave function. At 
this end, one can think about a suitable space which meets the properties of / and g; Lagrange 
interpolation functions are not suitable in this situation, since their first derivatives do not 
match the continuity property. So we consider instead a type of Hermitian functions (known 
as a generalization of the Lagrange functions) which are continuous and admit continuous first 
derivative. 

The boundary conditions need special concern, they play a crucial role of choosing the space 
of discretization; Since the wave functions are assumed to vanish at the boundaries and by the 
smooth property of these functions, the way they move toward the boundaries should be in 
damping manner, i.e with vanishing velocity, this implies zero derivative boundary conditions 
should be considered as well (except the case when k = ±1 at the lower boundary, see Remark 
1 below). Physically this is clearly reasonable, since the electron is neither expected to be close 
to the nucleus nor escaping to infinity. 

Cubic Hermite (CH) interpolation functions turn out to be sufficient to fulfill the requirements. 
Such functions are third-degree piecewise polynomials consisting of two control points and two 
control tangent points for the interpolation. That means there is a control for both the function 
values and the derivatives at each nodal point Xi. 

To study CH functions, let us first introduce the following spaces 



(i) For the states lsi/2 and 2pi/2 (k = — 1 and 1 respectively), the boundary conditions for 
the derivative of the components of the wave function are partially different, specifically 
at the lower boundary. Le if dQ^^ and denote respectively the upper and the 

lower boundaries, then v'Iqqup = and v'lg^io ^ 0. This is due to the fact that the 
corresponding wave function do not vanish in a damping way near the origin, see [20] for 



. J{o{n) = {v£:K{n):v'\9n = 0}. 



Remark 1. 
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more details. Thus, for k = ±1, the same functional space !Ko(r2) is considered but with 
small modification on the functions derivative at the lower boundary. Here we will keep 
the same notation !Kq{Q) for the space for all k's, but when we mean the cases k = ±1 
the right derivative condition at 90'° should be considered, 
(ii) For the sake of simplicity and as a matter of comparison, in the following computations 
of the energies of the electron in the Hydrogen atom, we do not use the right lower 
boundary conditions for k = ±1 as stated above. Instead we just assume zero for the 
derivative of the components of the wave function at 90'°, where the result might be 
slightly changed but does not affect the essence of the comparison. Also, without loss of 
generality, from now on we will assume v'\gQ = for all k. 

Let be the finite dimensional subspace of Jio on the partition kh spanned by CH basis 
functions. To summarize, possesses the following properties: 

(z) It is a set of continuous piecewise CH polynomials. 
(a) \/v S V^, v' exists and piecewise continuous. 

[in) yv € V^, v\qsupp{v) = v'\gsupp{v) = Oi where dsupp{v) is the boundaries of support of v. 
(iv) It is a finite dimensional vector space of dimension 2n with basis {4>j,i}^=i and {</'j,2}j=i 
given below. 

To approximate a function Uh G V^, where the same partition kh of the same distribution is 
considered as before, Uh can be written as 

n n 

(16) '^h = ^ ij4>j,i + X] , 

and are the unknown value of the function and its corresponding derivative at the nodal 
points Xj respectively, and and (j)j^2 are the basis functions of the space having the 
following properties 



1, Ifj = ^, 
, Otherwise . 



and 



'^i,2(^i) -jo! Otherwise. 



4>'j,i{xi) = (pj,2{xi) = Vi = 1,2, 



n. 



It follows from the conditions above that (/>j^i interpolates the function values whereas (f)j^2 is 
responsible of the function derivatives at the nodal point Xj. For non-uniform mesh, <j)j^i and 
(f)j^2 are given by the following formulae (see also Figure 1 below, where the two basis functions 
are depicted for uniform and nonuniform meshes) 



4-(x — Xj-i)"^ — -rjix — Xj_l)^(x — Xj) , X ^ Ij , 

j 

j;t^{x - Xjf + - Xjfix - Xj+i) , X G Ij+i , 



(IT) few = ^ ,^_^(,_ 

r ^{x - Xj^i)'^{x - Xj) , X £ Ij , 

(18) <Pj,2{x) = I (^_^.)__L_(2;-Xjf + ^(x-rEj)2(x-X,- + l), + 

The approximation error using CH basis functions in the subinterval Ij is given by 

(19) |m - < ci/i^||n('^)||2,oo(/^.) , 
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Figure 1. CH basis functions with uniform distributed nodal points(Left) and 
nonuniform distributed nodal points(Right). 



where ci = ggj, and h = maxj hj. So the error bound is obtained individually for each subinterval 
Ij, yielding a fine-grained error bound, which means that CH basis produces more accuracy 
compared to the linear or quadratic interpolation function in general. 

To construct FEM for the radial Dirac equation using CH basis functions, we as usual multiply 
([2]) and ([3]) by test function v € !Ko(il) and integrate over Q. To discretize the system we assume 
/ and g are elements of , thus they can be written as 



(20) 
(21) 



9(.x) = ^Ci</'i,i(a;) + X]^i'^J'2(x) , 
i=i i=i 

where Q and Cj sue the nodal value and the nodal derivative of / respectively at xj, and and 
are the corresponding ones associated to g. This yields 



CK, 
X 



+ 



n n 

^ {w^{x)(t)j^i{x) , v{x)^Cj + ^ {w^ {x)(t)j^2{x) , v{x)^C'j 



n n 



(23) 5^ (c(/<;. i(x) + ^-^4>,^i{x) , ^;(x))0 + Y ('^'^^■,2(^) + v'^:'-.2(x) , ^(x))cj + 



CK 
X 



n 



+ Y\^ (3;)</'i,i(a;) , v{x)j^j + Y\^ {x)4>j,2{x) , w(x)jCj- 
j=i i=i 

n n 

i=i i=i 



10 HASAN ALMANASREH, STEN SALOMONSON, AND NILS SVANSTEDT 

Let v{x) be an element of V^, and consider (j22p and (|23p first with v{x) = 4'i,i{x) ^^'^ then 



with v{x) 

(24) 

where 

(25) 

and 

(26) 



],2{x)- This yields the following system 

AX = X'BX , 



A 



mc^MMooo + MM^QQ 


-cMMqiq + CkMMqqi 


cMMoio + ckMMqqi 


-mc^MMooo + MM^QQ 



MMooo 








MMooo 



The vector X is the unknown given hy X = (CjC'iCjO) the general block matrices MM^^^ 
are defined as 



(27) 

where 
(28) 



MM^rst 





Ksta,2) 




^rst{2,2) 



Tables 3 and 4 contain the first six computed eigenvalues of the radial Dirac operator for the 
Hydrogen atom, with n = 100 and 400 interior nodal points using CH basis functions. The 
computation is run for k = — 1, 1, and the right column represents the exact solution for k = —1. 

Table 3. The first six computed eigenvalues for the electron in the Hydrogen 
atom using CH basis functions with 100 nodal points. 



Level 



1 



-1 



Rel. Form, k 



-1 



1 -0.50000632471 

2 -0.12500207951 

3 -0.05555629341 

4 -0.03125018386 

5 -0.01982545837 

6 -0.01085968925 



-0.50000665659 
-0.12500207951 
-0.05555629338 
-0.03125018404 
-0.01982545886 
-0.01085968695 



-0.50000665659 
-0.12500208018 
-0.05555629517 
-0.03125033803 
-0.02000018105 
-0.01388899674 



Table 4. The first six computed eigenvalues for the electron in the Hydrogen 
atom using CH basis functions with 400 nodal points. 



Level 



1 



K 



-1 



Rel. Form, k 



-1 



1 -0.50013790178 

2 -0.12500208021 

3 -0.05555629517 

4 -0.03125027925 

5 -0.01985891281 

6 -0.01116648473 



-0.50000665659 
-0.12500208018 
-0.05555629518 
-0.03125027916 
-0.01985888664 
-0.01116629119 



-0.50000665659 
-0.12500208018 
-0.05555629517 
-0.03125033803 
-0.02000018105 
-0.01388899674 
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It is noted, from the tables above, that some instilled spurious values are removed (the values 
that were present between level 3 and level 4). Also the speed of convergence to the exact 
eigenvalues is enhanced as the number of nodal points is increased. Unfortunately, part of the 
instilled spuriousity is still present for most values of Z, also the coincidence remains unsolved. 

The unphysical coincidence phenomenon assigns almost the same energies for 

N i(k = 1) and N i (k = -1) , N > 2. 
N 3(k = 2) and N 3(k = -2), N > 3. 
iV 5 (k = 3) and iV 5 (k = -3) , iV > 4. etc. 

The occurrence of this phenomenon is deeply studied for both nonrelativistic and relativistic 
cases; In [18], the coincidence of energies is proved for the same values of k that differ in sign 
via studying the commutation of Dirac operator with Biedenharn-Johnson-Lippmann (BJL) 
operator in the relativistic case. Also in the nonrelativistic case, the energy dependence on the 
quantum number ^ = n + |«;| is proved, which implies the energy independence of the sign of k. 
The coincidence of the energies in the finite space is also studied in [19] , where the spuriousity 
in general is interpreted as an effect of the same treatment of both components of the wave 
functions. 

As it is known that the exact solution of the Dirac operator with Coulomb potential for point 
nucleus results in different lowest bound energies for different values of k. In this work, as it is 
pointed before, we relate the problem of energies coincidence to the numerical implementation. 
Roughly speaking not to the method of approximation, but to the proposed spaces of discretiza- 
tion. 

In the last computations we imposed the test functions to live in the same space as well the 
trial functions, that is the usual Galerkin method. As we have seen, this results in a solution 
not cleaned from spurious values. However, it is well-known that the Galerkin method when it 
is applied to convection dominant problem, the solution will be upset by perturbations, which 
is getting worse with the increase in the convection size. 

Nevertheless, it is assumed non uniform mesh (exponentially distributed nodal points) to 
match desirable requirements of high resolution near the nucleus, where the wave functions 
oscillate rapidly compared to their oscillations in a region away from it. This means that for 
each nodal point Xj there are two adjacent systems of what are called fine- mesh grid with 
much larger coarse mesh. Hence when the wave function crosses the interface between these 
two regions, its phase is altered to fit the unbalanced change in the displacement size. One 
can understand the concept by regarding the variant mesh as different media to the generating 
waves, where most of those waves are not resolvable in two different meshes at the same time. 
We refer to [5] and [HI [12] for more details. 

Also, from numerical algebra point of view, one considers the linear system given by (j24p and 
posteriorly notes that the sign of k that appears as a factor of the block matrix TWMqoi does 
not contribute in determining the eigenvalues, which is entirely incorrect from physical point of 
view. So what is needed is to let the sign of k play a role in eigenvalues definition. This can be 
achieved by clever and justified addition of terms that includes k without deforming the original 
equations. These motivations suggest to use an alternative method to Galerkin formulation that 
does not demonstrate instability at the time treats the phenomenon of coincidence. 

Streamline Upwind Petrov-Galerkin (SUPG) method is used to solve the problem, which 
consistently introduces additional stability terms in the upwind direction, these terms are based 
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on the residual quantities for the governing equations and on the modification of the weighted 
test functional space. The latter is understood as adapting the test function t) from being {v, 0) 
and (0,u) to be {v,tv') and {tv',v) respectively, so it is a type of residual corrections added 
to the original equations. Tau ,t, is called the stability parameter which we are investigating, 
where its derivation is the main part of the upcoming section. 

2. Weighted Functional Space 

To stabilize FEM approximation applied to the Dirac operator, modified SUPG is used to 
formulate the problem. This consists in adding suitable stability terms to the standard Galerkin 
method. The SUPG method is designed to maintain the consistency, so that the solution of the 
original problem is still a solution of the variational equations. 

The idea behind SUPG is to introduce a diffusion term (n' , v') which eliminates the instability 
and enhances the approximation without modifying the problem. Several approaches can be 
implemented to create such term. To mention, we can just artificially add {au' , v'), where a 
is a constant that controls the diffusivity size, this method is first order accurate at most. Or 
the artificial diffusion term can be added in the direction of the streamlines to avoid excess 
diffusivity [H |4j, even though this method introduces less crosswind diffusivity than the first 
mentioned, but it is still inconsistent modification. The methods mentioned above result in a 
modified problem differs from the original by the addition of the terms which alter the structure 
of the problem and force the exact solution to be no longer satisfying the variational equations. 

To formulate the method, consider the radial Dirac equation 

Hr(p = X(p , where tp = (/(x) , g{x)Y and 
H-f ^^(^) + 



which is equivalent to 
(29) 



w+{x)f{x) - cg'{x) + fg{x) \ _ fix) 
cfix) + ffix) + w-{x)g{x) J g{x) 

Define the residual functional of each equation as 

(30) Rl{f,g){x) = w+{x)f{x) - cg'ix) + ^gix) - A/(x) = {W+ f - eg' + ^g){x) , 

(31) Rl{f,g){x) = cfix) + ^fix) + w-ix)gix) - Xgix) = iW-g + cf + ^/)(x) . 
Here W^ix) = w^ix) - A. 

The previously derived Galerkin discretization with CH basis functions reads 

(32) / V^Hripdx = X / oVc^a;, 

Jn Jn 

where is (v,0)* and (0,t;)*, and 



(33) ^ix) 



m ] ^ ( E-=iO<^„i(^) + E;=iCj<^i,2(x) 

So far with Galerkin approximation the components of as well as / and g are elements of . 

SUPG is formulated based on modifying the test function d to a form that includes v' as a 
correction term to introduce the required diffusivity. Hence we assume v € as well / and g, 
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but d ^ i^h)"^ is just continuous function. I.e, let D be {v,tv'Y and {tv' ,vY in ([32]) . where r is 
the stabihty parameter to be studied soon. This leads to 



(34) 



[w+f ^v) + {-cg' +''-^g,v) + {Rl{f, g) , tv') = X{f , v) 



and 



(35) 



{cf + 2^f,v) + [w-g , v) + {Rlif,g) , rv') = X{g , v) 



Each of the discretizations above, using the new weighted test functions, is the usual Galerkin 
formulation with additional perturbation terms consist of the weak variational form of the resid- 
ual of the opposite equation with basis function tv' . This keeps / and g, the exact solution, 
satisfying the weak formulation without modifying the problem. 

In matrix notations, the system AX = X'BX is obtained as before, but A and 23 are slightly 
perturbed by additional matrices factored by r 



(36) 



A 



mc^MMooo + AfMo^o+ 
+ctMMiio + ctkMMioi 


-cMMqio + ckMMooi + 
-mcVMMioo + tMMYoq 


cMMqio + ckMMooi + 
mcVMMioo + tMMYoo 


-mc^MMooo + MM^Q^+ 
-crMMiio + ctkMMioi 



and 



(37) 



MMooo 


tMMioo 


tMMioo 


MMooo 



The unknown vector X and the generalized block matrices MM^^^ are as defined before. It is 
notable from the system above that the resultant block matrices A and S are not symmetric any 
more, in this situation complex eigenvalues may will begin to appear, which of course what we 
should avoid in the computations. To be more precise, the appearance of complex eigenvalues 
depends on the size of r, where they do appear for large size. For small size of r one can consider 
the above system as the usual system that corresponds to the Galerkin approximation (which 
is symmetric) with an addition of small perturbation of size r, which still admits real eigenvalues. 

Now, the main task is to determine the stability parameter r that completes the scheme of 
removing the spuriousity for both categories and improves the convergence. The derivation r 
assumes non full dependence on the exact solution of the complete operator for point nucleus, 
instead limit operator is assumed. Parallel with considering the dominant terms relative to 
the speed of light. Before proceeding into details, we will give some lemmas without complete 
proofs, where the proofs in some cases are simple. 

The following lemma provides the approximated values of the radial function / and g at the 
nodal point Xj, where backward and forward derivative approximations are implemented, hence 
the error is 0(/i). 
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Lemma 1. For the Dirac spinors f and g, let Cj-i? Cj+i; Cj-ij be the f 's and g's 

nodal values at xj-i and Xj+i respectively, then the following holds 

^ (l + + (^( - + V(xi)) - '^xy, ^ 

e,_, = (l - + ( - '^{mc' + V(x,fj + '^X)(, ^ 

Proof . Consider the two-equation system of the radial Dirac equation 
(mc^ + V{x))f{x) + c{- g'{x) + ^g{x)) = \f{x) 

and 

c{f'ix) + ^/(x)) + ( - mc" + V{x))g{x) = \g{x) . 

Assuming the above system for arbitrary Xj E kh, and using the backward and the forward 
difference approximations for the derivatives (backward =^ f'U — -^^-^^'^ /(^i-i) _ Cj-i ^^^^ 
forward =^ /'U — ■^^^■'+^^ -^^^^^ = ^i±i— ^) one gets the desired results. ■ 

For the computed matrices MMqqq, MMiqq, MMqio, and MMhq in the block systems p6p 
and (j37p . the exact element integrals are obtained by the following lemma. For the remaining 
matrices one can calculate the exact element integrals, but it is rather hard to get them simplified. 
Therefore, we just point out in Remark 3 notations for the desired values without writing the 
explicit expressions. 

Lemma 2. The following table is the exact element integrals for some matrices in the generalized 
system. 

Table 5. The exact element integrals for some matrices in the generalized system. 



-^.^^^^ Index 


-....(J olurnn 
Row ^^^■"^----.^^ 


j - 1 


j 


i + 1 


j - 1 + n 


j + n 


3 + l + n 


Affl/ooo 


j 








420 "j + 1 


210V'j+l "jj 


420 "j + 1 


j + n 


-120 " j + 1 


2To(''?+i -''}) 


-120 " j + 1 


140 "j + 1 




140 "j + 1 


MMioo 


j 

j + n 


1 
2 





1 
2 


lT)''j+l 
60 "j+1 


- ll)('»j + i + hj) 



TO^'j + l 
60 "j + 1 




Tlj(^j+i + hj) 


-To^j+1 


JWMoio 


j 


1 

2 





1 
2 


^ TO ''j+1 


j^(?ij + i + hj) 


- TO'lj + l 


j + n 


Ti7''j + 1 


-TTiC'j + i +''-;) 


TFT + 1 


60 "j + 1 





60 "j + 1 


Amino 


j 


6 1 
5 fej + i 


5 hj^ihj 


6 1 
5 hj^i 


1 

10 





1 

10 


j + n 


1 

10 





1 

10 


- i^''j + l 


^{hj + i +hj) 


^ro''j+i 



Proof . The proof is straight forward by evaluating the integrals. 
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Remark 2. The basis functions consist of two parts, one corresponds to the function value and 
the other to the function derivative (the latter with no considerable contribution to the function 
values) at the nodal points. Therefore we will, for simplicity, just take into account the part of 
the basis functions that contributes to the function values at the nodal point only. Thus, the 
upper left (shaded) three-cell corner of each matrix of the above table is considered. 

Remark 3. For the other matrices in the block matrix A, iVlMooi and M/Vf loi (where MMq^q 
and MVlYoo can be written respectively as —ZAIMqqi and —ZNIMiqi for V{x) = — ), we will 
use the following notations for the calculations of the element integral using the part of the 
basis functions that contributes only the function values at the nodal points as indicated in the 
remark above. Namely as a matter of notation we denote the following 

Table 6. The element integrals notations of the matrices MMqqi and iVlMioi 
for the j^^ row. 



Index 


j-1 




J + 1 


Matrix 


3 


MMooi 








MMioi 









Now we are at the position to state the main theorem of the stability parameter r. 

Theorem 1. The mesh- dependence stability parameter r that appears in the formulations 
and ( 135|) is of the following form 

(38) T .- Tj = —hj+i 



35 '^\hj+i + hj) 

Before proceeding, we introduce the following notations to ease handling the proof. 

{hj+i + hj) 



ci 



2c 



9 hh Z Z 

C3 = ——hj+i{hj+i - hj) - KSj-i{hj+i - hj) - ^^(^j+i + hj)Tj + -{rj+ihj+i + 

-rj~ihj)Tj . 

6c(hj+i-hj) m'^c^ , , 1.^^ 2n/ , , n 

" y hj+^hj ^ (^i+i + ^■) + — ( V ~ ^'^ - ^i-i^i) • 

C5 = — ^—{hj+i + hj) - mcZ{rj+ihj+i - rj_ihj) + 5 ~ ~ ^3) + CK{rj_i + 

+rj + rj+i) . 
C6 = mc^ — ihj+i - hj) . 

Z 9 TTtC^ Kj 

C7 = -Z{2sj-i + Sj) + mc^KSj_i(/ij+i - hj) + - — (/ij+i + hj) - — /ij+i(/ij+i - hj) . 

9 

C8 = -—hj+i{hj^i-hj) . 
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t\j Z s 9 Z 

cio = ^(^+1 + ^J>i ^(^i+i - hj) + K(rj+i/ij+i - rj-i/ij)Tj - — — hj+i{hj+i + 

6mc2 1 
6 1 

Ci2 = -T — 7 - ^j) + mc^K(rj+i/ij+i - rj^ihj) - Z{rj_i + + rj+i) + 



+ ——{hj+i + hj) . 



2x, 

ci3 = 7^"^^c^^j+i(^i+i - ^i) H —{hj+i - hj) . 

9 mcZ CK , , 

ci4 = T^^^^'^J+iv^i+i ~ "-j) + CK(2sj_i + Sj) - + ^i) • 

Z s Z 
ci5 = —{mc^ H - hj) . 

C Xj 

ci6 = -^^^(mc^ - — - hj) . 

C Xj 

The following lemma provides the behavior of the eigenvalues in the vicinity of x at infinity. 
Lemma 3. Define the operator 

y cDx —mc^ 
Then for the radial Coulomb-Dirac equation 

+ ' V{x) ))[ g{x) ; - ^ V 9i^) J ' 

the only accumulation point of the eigenvalues X is mc? . 

Proof . See [9]. ■ 

We now give the proof of the main theorem. 
Proof . Consider the weak formulations (134p and (135p . rewrite both of them as the following 



matrix-system 

(39) (mc^ - A)MVfoooC - cMMqioC + ckMMqoi^ - ZMMqqiC + ctMMuqQ + 

+CKrMMioiC - {mc^ + A)tMMiooC - ^rMMioiC = 

and 

(40) (mc^ - A)rMMiooC - crMMno^ + cktMMiqiC - ^rMMioiC + cMMqioC + 

+ckMMooiC - (mc^ + A)MMooo? - ^MMooi^ = . 
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Where C = (Ci,'" ,Cj-i,Cj,Cj+i, ■ ■ ■ ,Cn,C'i,--- , Cj-i, Cj, Cj+i, ■■ ■ ,C) and ^ = (Ci,-" 
Cj+i, ■ ■ ■ , Cn, Cii ■ ■ ■ ; Cj) Cj+i; ' ' ' i C'n) ■ To get T locally, that is Tj, for each subelement of the 
mesh, we consider the above equations for arbitrary j cell. Employing Remark 2 and Remark 3 
together with Lemma 2 we end up with 



(41) (mc^ - a) (^/ij+iO-i + ^(^i+i + hj)Cj + ^hj+iCj+i) - c(- ^j-i + + 



6 1 6 (/ij+i + /ij) 6 1 \ / \ 



and 



(42) Tj I mc~ 



^0-1 - ^0+1 



6^_^^^ + 6 (/tj+i + ^i) ^ 
5 /ij+i 5 hj^ihj 



6 1 



5h 



■j+i 



+CK(sj_iCi-i + SjCj + Sj-iCj+i) - (mc^ + a) + ^(^i+i + ^i)0+ 



0. 



Using Lemma 1 (to substitute the nodal values Cj-i) O+ii Cj-i) and Cj+i)) equations above 
are written as 



(43) (I + CKs,_i - ^(mc2 + A)r, - Zr^^ir,) ((l - + (- ^(mc^ " ^) + ^ 

^ /ij)("ic2 - A) - Zsj + y ^^^7^^;^^^^i + cK^i^i) (o) + (c/^Sj - ^r^Tj) + 
+ (- I + -.,-1 + i(mc^ + A)r, - Zr,+,r,) ((l + + {^{mc' - ^) - ^A)o) + 



;A/^,+,(mc^ - A) - - |:^r, + c.r,.,r,) ((l + + (^^-"^^^ 



-mc 



-|A)e,) + (4/..>i(-c^ - A) - - l^r, + c.r,^,r,) ((l - ^)0+ 



+ (-^(-mc2-|) + %iA)e,)=0 
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- + CKSj^i + -{mc^- X)Tj-Zrj^iTj)( (l + ^)Cj + (^( 



-mc 



-)--A)6,| + 



+ (- ^(hj+i + hj){mc^ + A) - Zs 



6c (/ij+i + hj) 



5 hj^ihj '■' 



+ (| + --.-1 - ^(-c^ - A)r, - Zr,H.ir,) ((l - + (- ^(-mc^ - |) + 



+ ^A)Cj ) + (- --/ij+i(mc2 + A) - Zsj-i + ^-r^Tj + CKrj^iTj ) ( (l 



6c 1 



70 



5 h 



+(-^( 



mc 



9 



6c 1 



+c.r,^ir,) ((1 + ^)^, + (^(mc^ - ^) - ^A)^) =0. 



Rewriting ()33|) and by collecting the terms of and of respectively gives 



(45) y\^—hj+i{mc-\)-Zsj^i 



6c 1 



5 h 



Tj- j (2 H {hj - /ij+i) ) + CK(rj_i + + rj+i + 



Kr,_i , Kr.i+i . 6c + /i,) 13,, , n/ 2 <\ / 2 

- ^/i^.+i)r,- + -^-J±—^T, + + h,){mc^ - A) - + (mc2+ 



|-A 



- KSj-ihj + -^(mc + A)r, -I ,6,,, 

2 2c c 



2 + KSj-ihj+i + ^— (mc + 



+A)r, 



i+i 



+ 



z 



Zs 



mc'^ H 1" A ) ( —hj+ihj(mc^ — A) H ^^^^-^/i,-!- 

/ \ 70c c 



2 

+ck(2sj_i + Sj) - Z(rj_i + rj + rj+i)r,- - - — {hj + /ij+i) H ^(/ij+i - /ij) + 



2xj 

^ / 2 \\/7 T\ ^^i — l^l 

+- — {mc + X){nj^i + hj)Tj -\ /i,r. 



_ Zrj+iK 

„ "■j+'i-'j 



6 = 



and 



(46) -Z{rj^i + rj+rj+i)Tj + CK{2sj-i + Sj) + -^(mc^-A)(/ij+i + hj)Tj - ^{hj^i + hj) + 



+ 



CK Sj-\ 



Zs 



{hj - hj+i) 
6 h 



Zrj_iK Zvj^iK 



Z 



/ 9 



+ — — /ij - ^j^rj - KVj.ihjTj - —hj+i{mc + A) 



Zso_i 6 , 

——hj+i + -Tj + Krj+ihj^iTj 



+ 



9. - 2 " 6c 1 \/_ K 



^-h,^,{mc^ + X)-Zs,.,+ ^^^^^ 

_ Qc {hj+i + hj) ^ 13 
^' 5 hj+^h, ^ 35 



Tj\\2-\ (/ij+i - hj)] + CK{rj_i + rj + rj+i+ 



^/., + ^/.,.i)r, - -4Sl!m^r, - ^(/.,+i + h,){mc' + A) - Z., + • • • 
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+ mc'' H h A 



hj hj 2 , \ 

— - nsj-ihj -^[mc - \)Tj + 



hjTj + 



'j+i 



+ 



V+1 
2c 



mc 



Zr 



6 = 0. 



Gathering the factors of A^, A, Tj, and the free terms in each equation for ("j and respectively, 
and using the defined above notations Cj's, one can simphfy ()45p and ()46p as follow 



(47) 



and 
(48) 



CiTjA^ + (C2 + C3)A + (C4 + C5)rj + (C6 + C7) 



c,+ 



CsA^ + (Cg + Cio)A + 



+ (cil + Ci2)Tj + (Ci3 + Ci4 + C15) 



6 = 



CsA^ + (cg - Cio)A + (-C11 + Ci2)rj + (-C13 + Ci4 + Cie) 



+ 



+ (C2 - C3)A + (-C4 + C5)Tj + (-C6 + C7) 



6 



0. 



We consider the case where major part of the difficulty of solving the radial Dirac operators 
comes in. The above formulation is reduced to the operator T given in Lemma 3, the limit 
equation at infinity. One can understand the issue as the derived Tj should guarantee the 
stability of the computations in the entire domain, particularly for large x, which is the operator 
T in one hand, and to consider the dominant part of the operator which causes the instability 
in the computations in the other. These motivations allow to consider (j47p and (j48p in the limit 
case 



(49) 



i + l + ^j)^,y2 



2c 



-(/.,+!- MA +(-^^—^ 



+ 



2 S 

m c 



+—mc^{hj+i - hj] 



+ 



9 



2 

6 1 



{hj+i + hj) jTj + 
{hj+i - hj)Tj\ + 



6mc 1 



9 



5 h 



■j+i 



{hj+i - hj)Tj + —m c hj+i{hj+i - hj) ^j = 



and 
(50) 
9 



— /ij+i(/ij+i - hj)\^ 



6 1,, , s , 6?Tic^ 

[hj+i - hj)Tj\ + 



m^c^hj^i{hj^i — hj) 



70 



+ 



+1 



5h 

-hA 



5 h 



-j+i 



-(/ij+i - hj)Tj + 



2c 



-r,A' 



9 

70 



{hj+i - hj)\ + 



6c (/ij+i — hj) 
5 hj^ihj 



+ 



1 

m c 



{hj+i + hj)j Tj - T^rnc^ {hj+i - hj )J ^j = . 

■ J — and p = —9/70. Divide 1)39]) and (j50|) by the quantity 

/ij+i — hj (t^ for non-uniform mesh). In the vicinity of c at infinity one gets the following 
dominant equations 



Let m = 1, and define Vj 



{hj+i+hj) 



(51) 
and 
(52) 



[pA - aj]Cj + [djX - bj]Cj = 
[djX + bj]Cj + [pA + aj]^j = , 
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where aj = -(f + t^j^j + pc^, hj 



6c^ 1 



TJJ - and dj 



5 h 



6 1 
5 h 



■To . 



— Qj djX — bj 
djX + bj pX + 



Equations (jSip and (j52p can be written as 
(53) 

Since ("j and are not identically zero for all j, then 
(54) 

which gives 
(55) 







pA — Oj djX — bj 
djX + bj pX + Oj 



0, 



Ai,: 



Since is the accumulation eigenvalue (Lemma 3, with m = 1) we will only consider the positive 
part of A above named as Ai. Now we would like to have |Ai — = 



lAi 







Ai Jl 36 1 



C [p 



25 hf^ 



1 ■' 



36c' 



1__,2 I c6,-72^2 I 6c4 



] 3 



5 hj+ihj J j 



12c3 



-hJ^pTj - pC VjTj - -25- h^Jj 



rf+ 



-rp c -r 'j 4900 ""i+i ' 



keeping in mind the c limit at infinity, the above formulation gives 
(56) -V^rf - -^/il^i = . 



4 J J 



4900 



The desired result is then obtained straight forward after substituting the value of Vj as defined 
before, and this ends the proof. ■ 



The derived r provides complete cleaning of spectrum pollution for both categories. Also it 
is notable that the expression of r treats the difficulty of the wave transferring between any two 
adjacent unbalanced mesh steps. The size of r is proportional to the mesh size, i.e since we are 
dealing with exponentially distributed nodal points, r has small size near the singularity x = 
due to the small mesh size, where it takes relatively large values in the region away from the 
origin which is dominated by coarse mesh. 

Tables 7, 8, and 9 show the first computed energies for the electron in the Hydrogen-like 
Magnesium ion for both point and extended nucleus with k = j2|. Table 7 shows the computed 
eigenvalues using the usual Galerkin formulation with linear basis functions. The number of 
interior nodal points used is 400. Table 8 shows the same computations using the stability 
scheme. Table 9 represents the computed energies for extended nucleus using uniformly dis- 
tributed charge with interior nodal points 397, where 16 nodal points are considered in the 
domain [0 , R] [R is the radius of the nucleus). 
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Table 7. The first computed eigenvalues for the electron in the Hydrogen-like 
Magnesium ion using usual FEM with linear basis functions for point nucleus. 



Level 


K = 2 


K = -2 


Rel. Form, k = —2 


1 


^18.0086349982 


-18.0086349982 


-18.0086349982 


2 


^.00511829944 


-8.00511829944 


-8.00511739963 


3 


-4.50270135222 


-4.50270135225 


-4.50269856638 


^ 


^2.88546212211_ 


-2.88546212205 


Spurious Eigenvalue 


4 


' -2.88155295096" 


-2.88155295095 


-2.88154739168 


5 


-2.00096852250 


-2.00096852249 


-2.00095939879 


6 


-1.47003410346 


-1.47003410350 


-1.47002066823 


^ 


Le1.13034880166 


-1.13034880167^ 


Spurious Eigenvalue 


7 


-1.1254569168f 


-1712545691683 


-1.12543844140 


8 


-.889228944495 


-.889228944484 


-.889204706429 


9 


-.720265553198 


-.720265553187 


-.720234829539 


^ 


Q.600492562625 


-.60049256262^ 


Spurious Eigenvalue 


10 


-.595258516248 


-.595258516277" 


-.595220579682 


11 


-.500185771976 


-.500185772005 


-.500139887884 


12 


-.426201311278 


-.426201311300 


-.426146735771 



Table 8. The first computed eigenvalues for the electron in the Hydrogen-like 
Magnesium ion using the stability scheme for point nucleus. 



Level 


K = 2 


K= -2 


Rel. Form, k = —2 


1 




-18.0086349985 


-18.0086349982 


2 


-8.00511739978 


-8.00511740020 


-8.00511739963 


3 


-4.50269856669 


-4.50269856719 


-4.50269856638 


4 


-2.88154739219 


-2.88154739270 


-2.88154739168 


5 


-2.00095939948 


-2.00095939991 


-2.00095939879 


6 


-1.47002066888 


-1.47002066924 


-1.47002066823 


7 


-1.12543844176 


-1.12543844201 


-1.12543844140 


8 


-.889204706068 


-.889204706109 


-.889204706429 


9 


-.720234827833 


-.720234827687 


-.720234829539 


10 


-.595220575840 


-.595220575531 


-.595220579682 


11 


-.500139880950 


-.500139880357 


-.500139887884 


12 


-.426146724530 


-.426146723650 


-.426146735771 


13 


-.367436809137 


-.367436807839 


-.367436826403 


14 


-.320073519367 


-.320073498169 


-.320073665658 


15 


-.281295132797 


-.281293164731 


-.281311119433 



To study the convergence property of the derived scheme, we compare the approximated 
eigenvalues of the electron in the Hydrogen-like Magnesium ion for point nucleus using the 
usual FEM as in Table 7, to those values obtained by the stability scheme as in Table 8. Ignor- 
ing the presence of the spurious values, one notes that the relative error in the approximation 
of the first 12 genuine eigenvalues using FEM is nearly 10~^. Whereas the relative error for the 
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Table 9. The first computed eigenvalues for the electron in the Hydrogen-like 
Magnesium ion using the stability scheme for extended nucleus. 



Level 


K = 2 


K = -2 


Rel. Form, k = —2 


1 




-18.0086349986 


-18.0086349982 


2 


-8.00511739975 


-8.00511740015 


-8.00511739963 


3 


-4.50269856673 


-4.50269856733 


-4.50269856638 


4 


-2.88154739230 


-2.88154739279 


-2.88154739168 


5 


-2.00095939956 


-2.00095940014 


-2.00095939879 


6 


-1.47002066903 


-1.47002066934 


-1.47002066823 


7 


-1.12543844179 


-1.12543844207 


-1.12543844140 


8 


-.889204706021 


-.889204706003 


-.889204706429 


9 


-.720234827640 


-.720234827433 


-.720234829539 


10 


-.595220575309 


-.595220574883 


-.595220579682 


11 


-.500139879906 


-.500139879215 


-.500139887884 


12 


-.426146722827 


-.426146721812 


-.426146735771 


13 


-.367436806543 


-.367436805088 


-.367436826403 


14 


-.320073514034 


-.320073492344 


-.320073665658 


15 


-.281294966822 


-.281292979627 


-.281311119433 



same group of eigenvalues using the stability scheme is not exceeding 3 * 10 . Thus, the speed 
of convergence is also enhanced. 

In Table 10, we provide the approximated eigenvalues for the electron in the Hydrogen-like 
Uranium ion using the stability scheme. The computations are obtained for different values of 
the quantum number k for extended nucleus. The number of nodal points used is 203 (13 out 
of them are used to discretize the segment [0 , B\). 

Conclusion. 

Our computations indicate that the SUPG scheme applied to solve the radial Dirac eigenvalue 
problem is stable in the sense of complete elimination of spectrum pollution. This approach is 
mainly compiled of two strategies; the first is the suitable choice of the trial functional space. 
The second is based on varying the test function to live in another space different from that for 
the trial function, this strongly depends on the derived stability parameter r. The derived r 
is a considerable achievement where its formula is rather easy to implement, and it yields full 
treatment of the spuriousity for both categories. 
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Table 10. The first computed eigenvalues for the electron in the Hydrogen-like 
Uranium ion for different energy levels using the stability scheme for extended 
nucleus. 



Level 




K = 1 


K = -2 


K = 2 


K = —3 


1 

X 












9 


1 9'i'i 0^89791 fi 


1 9'i7 997'^SR41 








Q 
O 


-ooo.uuioouyuo 


'i'^Q 0'^'^QQ0^9R 
-ooy.uooyyuozu 


1 OSQ fil 1 41 •i'i9 






A 


9QPi 078798090 


9Qf^ 9*^90/14 ^07 
-zyo.zozu^^ou ( 


/18Q n'^7n8Pil "^4 
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-4oy .UO ( uo4yOU 




o 


1 8^1 '^Q^OQOfi'^fi 
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-3.94434475273 


-3.94459787924 


-4.15983518735 


-4.15985667790 


-4.42238191278 


34 


-3.71370352567 


-3.71393684674 


-3.91041067459 


-3.91043403573 


-4.14944192144 


35 


-3.50263193199 


-3.50284782636 


-3.68268594092 


-3.68271127049 


-3.90093812765 


36 


-3.30896641306 


-3.30916694458 


-3.47420485278 


-3.47423225164 


-3.67402878591 


37 


-3.13083310725 


-3.13102007468 


-3.28284801467 


-3.28287758577 


-3.46627240831 


38 


-2.96660233349 


-2.96677731307 


-3.10677874375 


-3.10681058942 


-3.27556186817 


39 


-2.81485102305 


-2.81501541978 


-2.94439883084 


-2.94443304812 


-3.10007081630 


40 


-2.67433187309 


-2.67448701724 


-2.79431185210 


-2.79434854846 


-2.93820980064 


41 


-2.54395104802 


-2.54409831318 


-2.65529221336 


-2.65533164690 


-2.78858986482 


42 


-2.42276031374 


-2.42290074014 


-2.52626102582 


-2.52630405138 


-2.64999168906 


43 


-2.30995095161 


-2.31008292285 


-2.40627974326 


-2.40632804795 


-2.52134098372 


44 


-2.20475351060 


-2.20486906775 


-2.29457262587 


-2.29462526625 


-2.40170120650 


45 


-2.10615086328 


-2.10624270502 


-2.19049573109 


-2.19053698126 


-2.29029930213 
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